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Preface 


Tills  report  is  the  fifth  concerned  with  research  accomplished  in 
connection  with  Navy  Contract  Nonr-433(00) , between  Dunham  Laboratory, 
Yale  University,  and  the  Office  of  Naval  Research,  Department  of  the 
Navy.  In  this  report  is  given  a discussion  of  the  solution  for  a 
nonlinear  differential-difference  equation.  Methods  of  attacking 
equations  of  this  general  sort  are  but  imperfectly  developed,  and 
only  approximate  solutions  can  be  obtained.  The  particular  equation 
considered  here  may  apply  to  several  physical  phenomena  of  interest, 
and  the  mathematical  analysis  is  of  interest  in  itself* 

The  research  was  carried  on  by  W,  J,  Cunningham,  with  the 
assistance  of  J.  G.  Skalnikj  the  report  was  written  by  the  under- 
signed. 


W,  J.  Cunningham 
New  Haven,  April,  1954 
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Abstract 

Certain  physical  phenomena  appear  to  be  described  by  the  nonlinear 
differential-difference  equation 

dx(t)/dt  = £a  - b x(t  - x)j  x(t) 
where  a,  b,  and  t are  positive  real  constants.  This  is  an 
equation  of  growth  in  which  the  growth  rate  of  a quantity  depends  in 
part  upon  the  value  of  that  quantity  at  some  earlier  time.  Methods 
of  obtaining  exact  solutions  for  this  sort  of  equation  are  unknown. 
Approximate  solutions  for  the  equation  are  obtained  by  several 
analytical  methods.  Variable  x can  never  go  through  the  value  zero, 
and  thus  reverse  its  algebraic  sign.  If  x is  positive  and  the 
product  ax  is  small,  solutions  approach  a limiting  value  a/b, 
either  monotonically  or  with  a decaying  oscillation.  If  product  ax 
is  larger,  a steady-state  oscillation  of  definite  amplitude  occurs. 

If  x is  negative,  solutions  run  off  to  negative  infinity. 

Examples  of  solutions  for  particular  values  of  the  parameters 
are  obtained  with  an  analog  computer. 
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I*  Phenomena  and  Equation  under  Consideration 

There  are  certain  liatural  phenomena  in  which  the  magnitude  of 

some  quantity  increases  at  a rate  proportional  to  the  magnitude  itself. 

A simple  example  is  the  growth  in  population  of  an  organism,'*'  where 

the  number  of  new  individuals  appearing  within  any  given  short  time 

interval  depends  upon  the  number  of  individuals  present  at  the  beginning 
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of  the  interval.  Certain  chemical  or  neuclear  reactions  may  operate 
in  a similar  way,  with  the  rate  of  reaction  proportional  to  the  amount 
of  end  product  that  is  present. 

Phenomena  such  as  these  can  be  described  mathematically  try  the 
differential  equation 

dx/dt  = ax  (l) 

where  x is  a measure  of  the  product  in  question,  t is  time,  and  a 
is  a positive  real  constant,  the  relative  growth  rate.  This  equation 
is  sometimes  referred  to  as  the  equation  of  growth.  Its  solution  can 
be  written  as 

x ~ x exp(at)  (2) 

where  xq  is  the  value  of  x existing  at  zero  time.  Curves  repre- 
senting x as  a function  of  t have  the  well-known  exponential  shape 
of  Fig,  1,  where  a family  of  curves  is  shown  for  constant  xq  but 
several  values  of  a.  As  t increases  without  limit,  so  also  does  x 
increase  without  limit,  the  rate  of  increase  becoming  proportionally 
larger. 

There  are  certain  phenomena  which  appear  to  be  governed  initially 
by  an  equation  such  as  Eq,  (1),  However,  as  the  quantity  represented 
by  x increases,  some  effect  comes  into  play  which  reduces  its  rate  of 

1,  H,  Margenau  and  G,  II.  Murphy,  Mathematics  of  Physics  and  Chemistry, 
(Van  Nostrand,  New  York,  1943),  p,  33 

2,  W,  Jc=t,  Explosion  and  Combustion  Processes  in  Gases,  (McGraw-Hill, 
New  York,  1946),  p,  2G2 
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increase.  Ultimately  x may  approach  a constant  value.  This  might 

be  the  case  for  a colony  of  biological  organisms  in  circumstances  where 

the  food  supply  or  living  space  is  limited.  As  the  population  of  the 

organisms  approaches  the  maximum  which  can  be  supported  continuously, 

the  rate  of  growth  must  be  reduced  in  some  way.  One  hypothesis  that  has 
3 

been  suggested  is  that  the  members  of  tha  population  somehow  recognize 
that  crowding  is  soon  to  occur.  Upon  making  this  discovery  . they  try 
consciously  to  reduce  their  reproduction  rate.  Because  of  a finite 
gestation  time,  the  actual  birth  rate  is  not  lowered  until  a definite 
time  after  an  effort  has  been  made  to  lower  it.  An  equation  describing 
this  kind  of  operation  is 

dx(t)/dt  **  jja  - b x(t  **  t)J  x(t)  (3) 

where  a and  b are  positive  real  constants,  t is  a constant  delay 
time,  and  x(t)  and  x(t  — x)  are  values  of  x at  the  instants  t 
and  (t  - x),  respectively#  The  effective  growth  rate  in  Eq.  (3)  is 
|a  - b x(t  - x)j  , It  is  less  than  parameter  a by  an  amount  propor- 
tional to  the  value  of  x existing  at  the  time  (t  - x),  earlier 
than  the  time  t at  which  the  derivative  dx(t)/dt  is  evaluated* 
Another  example  where  Eq.  (3)  might  occur  is  in  the  control  ff 
some  reaction  which  fundamentally  is  governed  by  Eq.  (l) . In  order  to 
prevent  the  reaction  running  away,  with  catastrophic  results  perhaps, 
some  modification  is  intentionally  introduced  into  the  system  to 
reduce  the  reaction  rate.  The  controlling  mechanism,  which  senses 
the  rate  of  reaction  and  takes  steps  to  change  it,  requires  a finite 
time  to  operate.  If  the  delay  time  is  of  fixed  value,  x,  the  equation 
applying  to  the  system  is  Eq.  (3). 

3.  This  equation  appears  to  have  been  introduced  in  the  paper, 

G,  E.  Hutchinson,  ^Circular  Causal  Systems  in  Ecology”,  Annals  of  the 
New  York  Academy  of  Sciences,  £0,  221,  (1948). 
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Still  another  example  where  this  equation  might  arise  is  in  the 
determination  of  the  path  of  a controlled  projectile.  A mass  particle, 
falling  freely  in  a constant  gravitational  field  with  no  retarding 
effects,  follows  a parabolic  path.  Furthermore,  the  curve  of  height 
as  a function  of  time  also  is  parabolic,  since  the  horizontal  component 
of  velocity  remains  constant,  and  the  horizontal  position  is  propor-' 
tional  to  time.  An  example  of  the  relation  between  height  and  time  is 
shown  in  Fig.  2.  Coordinates  x and  t are  chosen  so  that  the  vertex 
of  the  parabola  is  located  at  x = 0 and  t = -t^,  For  convenience, 
positive  x is  plotted  below  the  t-axis«  At  zero  time  the  particle 
is  at  the  point  x * xq  and  t = 0,  The  parabolic  curve  can  be 
approximated  over  a part  of  its  length  by  an  exponential  curve 
suitably  chosen.  Such  a curve  is  shown  also  in  Fig.  2.  An  actual 
mass  particle,  falling  in  the  gravitational  field  of  the  earth,  is 
retarded  by  the  effect  of  resistance  with  the  atmosphere.  This  retal- 
iation ultimately  makes  the  vertical  component  of  velocity  constant, 
so  that  the  curve  of  Fig,  2 would  approach  asymptotically  a straight 
line  of  constant  slope.  The  curve  then  could  not  be  approximated  by 
an  exponential  curve. 

Instead  of  the  particle  falling  freely,  it  might  be  subject  to 
control,  the  intent  of  which  is  to  make  the  path  become  horizontal  at 
some  definite  value  of  x.  Again  the  control  system  requires  a finite 
time  to  operate.  If  this  delay  time  is  constant,  and  if  the  free  fall 
of  the  projectile  is  assumed  to  be  essentially  exponential  in  shape, 

Eq.  (3)  may  describe  the  path  which  the  projectile  follows. 

Equations  similar  to  Eq,  (3)  occur  in  economic  studies^  of 
business  cycles,  where  time  delays  occur  in  various  steps  of  the 
business  operations, 

4.  R.  M,  Goodwin,  Econometrica,  19,  1,  (1951) 


5 


In  general,  any  phenomenon  in  which  some  quantity  tends  to  grow 
at  an  ever-increasing  rate,  but  is  subjected  to  a throttling  effect 
involving  a constant  time  delay,  may  be  described  by  Eq„  (3) . 

The  analysis  cf  nonlinear  differential-difference  equations^  seems 
not  to  have  been  studied  in  great  detail,  so  that  methods  of  finding 
exact  solutions  for  Eq«  (3)  are  unknown.  Qualitatively,  however,  it  is 
apparent  that  if  x(t)  ever  vanishes,  the  value  of  dx(t)/dt  also 
must  vanish.  If  this  occurs,  x can  never  become  different  from  zero 
again.  Thus,  if  x(t)  is  not  zero,  it  must  always  retain  the  same 
algebraic  sign,  and  the  signs  of  both  x(t)  and  x(t  — x)  must  be 
identical. 

If  x(t)  is  positive,  the  sign  of  dx(t)/dt  may  be  either 
positive  or  negative,  depending  upon  the  relative  magnitudes  of  the 
terms  on  the  right  side  of  Eq,  (3),  Thus,  the  appearance  of  oscil- 
lations with  x(t)  positive  is  allowed,  A steady  value,  x(t)  = 
x(t  - x)  «=  a/b,  for  which  dx(t)/dt  * 0,  also  is  a possible  solution. 

If  x(t)  is  negative,  the  sign  of  dx(t.)/dt  always  is  negative,  and 
the  solution  can  only  go  to  negative  infinity, 

II,  Degenerate  Case  for  x = 0 

A simple  case  of  Eq,  (3),  and  a profitable  one  to  take  as  a 
starting  point,  is  that  for  which  the  delay  time  x is  zero,  so  that 
the  equation  is^ 

dx/dt  = a?.  - bx^.  (4) 

5,  R.  Bellman  and  J,  K,  Danskin,  "Stability  Theory  of  Differential- 
Difference  Equations,"  Proceedings  of  Symposium  on  Nonlinear  Circuit 
Analysis,  (Polytechnic  Institute  of  Brooklyn,  New  York,  1953).  p»  107. 
This  reference  contains  a long  bibliography, 

6,  This  is  the  Verhulst-Peari  equation.  See  A,  J»  Lotka,  Elements  of 
Physical  Biology,  (Williams  and  Wilkins,  Baltimore,  1923),  p,  64. 
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This  equation  can  be  studied  initially  by  considering  a plot  of  the 

new  variable  y s dx/dt  » x as  a function  of  Such  a plot  is  the 
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Poir.car8'  phase  plane,  widely  used  in  studying  second-order  equations* 
The  phase-plane  curve  for  Eq»  (4)  is  shown  in  Fig*  3* 

There  are  two  points,  x = 0 and  x = a/b,  where  y = 0*  These 
are  points  of  equilibrium,  since  no  change  can  occur  if  x initially 
has  either  of  these  values.  If  x initially  has  some  other  value, 
changes  will  occur.  If  y is  positive,  x changes  to  become  more 
positive,  and  vice  versa.  Thus,  changes  in  x with  respect  to  time 
take  place  in  the  directions  indicated  in  Fig*  3*  The  point  x = 0 
is  unstable  in  the  sense  that  x tends  to  run  away  from  this  point 
with  time.  The  point  x = a/b  is  stable,  since  x tends  to  converge 
toward  this  point.  Both  points  are  termed  nodal  points,  with  the 
curve  of  y as  a function  of  x approaching  from  a definite  direction* 
An  exact  solutior  for  Eq,  (4)  can  be  found,  considering  it  as  an 
example  of  Bernoulli’s  equation.  The  substitution  w = l/x  is  made, 
giving  the  linear  equation 

dw/dt  + aw  = b (5) 

which  has  the  solution 

w = b/a  + C exp(-at)  (6) 

with  C an  arbitrary  constant.  Solution  of  Eq,  (4)  is,  then, 

x = £b/a  + (I/xo  ~ b/a)  exp(-at)j"’'!'  (7) 

where  x = x at  t = 0, 
o 

The  nature  of  this  solution  can  readily  be  compared  with  Fig.  3* 

As  t becomes  infinitely  large,  x always  approaches  the  value 
x = a/b  monotonically.  If  xq  is  positive,  the  curve  for  x as  a 


7.  N.  Minorsky,  Nonlinear  Mechanics,  (J.  V/,  Edwards,  Ann  Arbor,  1947) 
Part  I 
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function  of  t is  continuous.  Its  curvature  is  zero  at  the  value  of 
t making  d x/d t **  0,  which  is  the  same  as  dy/dx  = 0,  This  value  is 
t = a ^ ln(a/hxo  - 1)  at  which  instant  x * a/2b,  and  the  slope  of 
the  curve  has  its  maximum  positive  value  at  this  point.  If  xq  is 
negative,  the  curve  for  x goes  to  infinity  for  the  value  of  t 
causing  the  quantity  in  the  bracket  of  Eq.  (7)  to  vanish.  This  value 
is  t = a ^ ln(l  - a/txo), 

A family  of  solution  curves  for  Eq,  (4)  is  shown  in  Fig.  4,  For 
small  values  of  x,  the  rate  of  growth  is  determined  primarily  by 
parameter  aj  for  large  values  of  x it  is  determined  primarily  by 
parameter  b.  The  ultimate  value  of  x depends  upon  the  ratio  a/b. 
These  solution  curves  apply  to  Eq,  (3)  with  t * 0,  and  with  both  a 
and  b positive.  Extension  to  the  case  of  reversed  sign  for  these 
parameters  is  self-evident. 

Ill,  Differential  Equation  Approximately  Equivalent 
III.  1.  Derivation  from  differential-difference  equation 

It  is  difficult  to  study  the  differential-difference  equation, 

Eq.  (3),  because  of  the  term  x(t  - x)  which  is  evaluated  at  a time 
different  from  the  other  terms,  x(t)  and  dx(t)/dt.  A differential 
equation,  with  all  terms  evaluated  at  the  same  instant,  is  more  easily 
analyzed.  The  term  x(t  - t)  can  be  expressed  by  the  Taylor’s  series 
x(t  - t)  « x(t)  - xdx(t)/dt  + (x^/2)dSc(t)/dt^ 

- (x^/6)d?x(t)/dt"*  (8) 

where  each  term  on  the  right  side  is  evaluated  at  the  same  instant,  t. 
If  the  delay  time,  x,  is  sufficiently  small,  its  higher  powers  become 
still  smaller  and  the  series  of  Eq.  (8)  may  converge  fairly  rapidly. 
Then  a good  approximation  to  the  value  of  x(t  - x)  may  be  had  from 
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only  the  first  few  terms  of  the  series.  If  either  x or  the  higher- 
order  derivatives  are  not  sufficiently  small,  the  series  converges 
slowly,  and  many  terms  must  be  retained  for  accuracy* 

If  only  the  first  three  terms  of  Eq*  (8)  are  substituted,  into 
Eq.  (5),  the  result  is 

x - a xx  + a x/bx  + a x = a p (.9; 

2 2 

where  a s 2/x  , p = a/b,  and  all  terms  are  evaluated  at  the  same 
instant.  If  x = 0,  Eq.  (9)  reduces  to  Eq.  (4)  as  it  should. 

Equation  (9)  is  a second-order,  nonlinear  differential  equation, 
representing  a system  with  a single  degree  of  freedom,  With  certain 
choices  of  parameters  its  solutions  are  oscillatory,  but  only  a single 
mode  of  oscillation  at  a single  frequency  can  occur  at  a given  time. 

A differential  equation  of  infinite  order  would  have  resulted  if  all 
terms  of  Eq.  (8)  had  been  used.  Such  an  equation  would  represent  a 
system  with  an  infinite  number  of  degrees  of  freedom,  having  the 
possibility  of  simultaneous  oscillation  at  an  infinity  of  frequencies. 
The  modification  of  the  differential-difference  form,  Eq.  (3),  to  the 
pure  differential  form,  Eq.  (9)>  has  brought  a simplification,  but  at 
the  expense  of  losing  the  possibility  of  simultaneous  oscillation  at 
several  frequencies.  Techniques  for  studying  a second-order  equation, 
such  as  Eq,  (9),  are  well  developed;  equations  of  higher  order  are  more 
difficult  to  analyze.  It  is  not  unreasonable  to  hope  that  solutions 
for  Eq,  (9)  will  be  similar  to  those  of  Eq,  (3),  although  exact  cor- 
respondence cannot  be  expected. 

Equation  (9)  itself  can  be  interpreted  as  applying  to  a physical 

system  somewhat  different  from  the  examples  used  in  connection  with 

2* 

Eq.  (3).  If  the  nonlinear  term,  a x/bx.,  of  Eq.  (9)  is  neglected,  the 


10. 

remaining  terms  are  linear,  and  the  equation  has  the  solution 

x = 3 + C exp(t/x)  cos(t/T  + 0)  (10) 

where  C and  0 are  arbitrary  constants.  This  solution  is  the  sum 
of  a constant  term  and  an  oscillation  with  amplitude  growing  with  time. 
The  nonlinear  term  of  Eq.  (9)  represents  damping  which  varies  inversely 
with  the  value  of  x.  Thus  the  solution  of  nonlinear  Eq,  (9)  might 
be  expected  to  resemble  Eq,  (10),  but  with  large  damping  occurring  at 
those  instants  when  x is  small.  It  is  not  unlikely  that  the  combi- 
nation of  this  damping  with  the  growing  exponential  factor  of  Eq.  (10) 
will  result  in  the  appearance  of  a limit  cycle, ^ an  oscillation  of 
fixed  amplitude  determined  solely  by  the  parameters  of  the  equation  and 
independent  of  initial  conditions.  Thus,  Eq,  (9)  might  apply  to  an 
oscillator  with  a linear  negative-damping  term,  and  a positive-damping 
term  varying  inversely  with  the  value  of  x.  Since  the  positive 
damping  would  become  infinite  if  x ever  goes  to  zero,  it  is  evident 
that  a steady  oscillation  can  occur  only  about  some  non-zero  mean  value 
of  x,  and  this  i3  provided  by  the  constant  term  of  Eq.  (9)» 

Equation  (9)  may  be  compared  with  the  well-known  van  der  Pol 
equation, 7 x - e(l  - x )x  + x = 0,  which  also  has  solutions  in  the 
form  of  limit  cycles.  In  the  van  der  Pol  equation,  positive  damping  for 
large  amplitudes  of  oscillation  occurs  symmetrically  at  either  extreme 
of  the  cycle.  In  Eq,  (9),  positive  damping  occurs  asymmetrically, 
being  large  only  for  instantaneous  values  of  x on  the  negative  side 
of  the  mean  value,  3* 


8. 
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N.  Minorsky,  ref,  7,  p«  62 
N.  Minorsky,  ref,  7,  p.  113 


ill,  2.  Singularity  at  x » g,  y = 0. 

The  qualitative  nature  of  the  solutions  for  Eq,  (9)  can  be  deter- 
mined by  studying  the  singular  points^  of  the  first-order  equation 
resulting  from  the  substitution,  y = x,  so  that  x » y = y dy/dx. 
This  equation  is 

dy  = a2L(T  - l/bx)y  - x + 3.1  /-jjn 

dx  y . K ' 

Its  solution  can  be  represented  as  a curve  on  the  phase  plane  having 

as  axes  the  coordinates  x and  y.  Singular  points  are  located  at 

those  values  of  x and  y which  make  both  numerator  and  denominator 

of  Eq.  (11)  vanish  simultaneously.  Only  one  singular  point  is  so 

determined:  its  coordinates  are  x = B and  y ■ 0.*  The  nature  of 

S 3 

the  solutions  near  this  singularity  can  be  found  by  replacing  x and 
y with 


x ^ x +u  = B + u 

s 

y — -*  yg  + v = v 

where  u and  v are  small  increments.  If  these  substitutions  are 
made  in  Eq.  (ll),  and  only  linear  terms  are  retained,  the  result  is 

dv  _ g2H(T  - 1/ a)v  - u] 

du  v * ^ ' 

Since  dv/du  = (dv/dt)/(du/dt)  = (l/v) (d2u/dt2),  Eq,  (12)  is  equivalent 
to  the  equation 

d2u/dt2  - a2(x  - l/a)  du/dt  + a2u  = 0.  (13) 

Solutions  for  this  linear  second-order  equation  depend  upon  the  roots, 

4-^  and  A.£,  for  the  characteristic  equation 

s2  - <x2(t  - l/a)  s + a2  = 0 . (14) 

10.  N,  Minorsky,  ref.  7,  Ch.  Ill 

* It  is  worth  remarking  that  a second  singularity  is  located  at  the 
origin,  but  that  this  point  is  not  a simple  singularity  and  cannot  be 
discussed  in  the  same  manner  as  the  first  singular  point.  It  is  con- 
sidered briefly  in  Sec.  III.  5# 


These  roots  are 


\p  \2  “ (ax2)  1 |(ax  - 1)  ± [(ax  - l)2  - 2a2x2J1/  2 j . (15) 

Solutions  are  norr-oscillatory  if  the  roots  are  real  quantities,  that 
is,  if  j(ax  - l)2  - 2a2x2J  > 0.  This  will  be  the  case  if 
(-2^2  - 1)  c (ax)  c (2^2  “1);  otherwise  oscillations  will  occur. 
Solutions  are  stable,  in  the  sense  that  they  do  not  increase  without 
bound  as  time  increases,  if  (ax  ) (ax  - l) < 0.  Actually,  of  course, 
parameters  a and  x were  assumed  positive  at  the  beginning  of  the 
discussion,  but  they  might  be  allowed  to  become  negative  under  seme 
conditions.  The  qualitative  nature  of  solutions  is  depicted  graphically 
in  Fig.  5.  The  notation  used  here  is  that  conventional  in  referring 
to  the  nature  of  a singularity,  A focus  refers  to  an  oscillating 
solution  with  amplitude  either  decreasing  or  increasing;  a node  refers 
to  a solution  approaching  a limiting  value  monotonically.  The  solution 
found  in  Sec.  II,  for  a>0  and  x = 0,  approaches  the  singular 
point  monotonically,  and  the  point  is  a stable  node. 

According  to  the  theory  of  Liapounoff,  the  nature  of  solutions 
for  the  nonlinear  equation,  Eq,  (9),  are  similar  to  those  of  the  linear 
equation,  Eq,  (13),  so  long  as  the  variables  u and  v are  suffi- 
ciently small.  Thus,  in  the  nieghborhood  of  the  singular  point,  solu- 
tions of  Eq,  (9)  have  the  properties  illustrated  in  Fig,  5« 

III,  3.  Approximate  solution  by  variation  of  parameters,  about  (x  ,y  ) 

S 8 

Still  more  information  about  the  solutions  of  Eq,  (9)  can  be 

12 

found  by  applying  the  method  of  variation  of  parameters.  This  method 
is  useful  where  an  oscillating  solution  occurs  with  only  small  changes 
in  either  amplitude  or  phase  taking  place  within  a cycle.  In  applying 

11.  N.  Minorsky,  ref.  7,  p«  51 

12.  N,  Minorsky,  ref.  7,  Ch.  X 


14. 

the  method,  it  is  convenient  to  rewrite  Eq.  (9)  as  a pair  of  first- 


order  equations. 


o o o 

y “ a (t  - l/bx)y  - ax  + a^p  . 


(16) 


(17) 


(18) 


If  the  term  in  y,  with  its  nonlinearity,  is  emitted  .from  the  second 
of  these  equations,  a generating  solution  is  found  as 
x = p + A sin(at  + 6)  - 3 + A sin  / 
y 3 aA  cos  / 

where  A is  the  amplitude  and  8 is  the  phase  angle  found  from  initial 
conditions,  and  / s (at  +8).  This  generating  solution  is  then  used 
in  the  complete  form  of  Eq,  (16),  allowing  both  A and  8 to  become 
time-dependent  instead  of  being  mere  constants*  The  result  of  this 
substitution  is 

A sin  / + A8  cos  / ~ 0 

O 

! j . a • j 3 , j or  A cos  / 

to  cos  - Ace  sin  t - a tA  cos  i - bgiTVTA/|3)  sin  ^ 

Solution  for  A gives 

A » a2A  cos2/  jjr  - a ^(1  + A/p  sin  4)  ^j» 

If  A/p  « 1,  approximately 

(1  + A/p  sin  /)  ^ » 1 - A/p  sin  / + A2/p2  sin2  / + • • • 

« 

and  A becomes 

A = a2vA  cos2  / - a2A/a  cos2  / + a2A2/aP  cos2  / sin  / 

- a2A^/ap2  cos2  / sin2  /+•••• 

If  the  amplitude  changes  but  slowly,  only  the  average  rate  of  change 
over  a cycle  need  be  considered.  Since 
cos2/  =*  ^(1  + cos  2/) 

cos2/  sin  / = T(sin  / + sin  3/) 

4 

cos2/  sin2/  = ^(l  - cos  4/) 


(19) 


15. 


« 

the  average  value  of  A over  a cycle  is 

(A)  « (ax2)-1  (ax  - 1)  A - b2A3/4a3x2 

»V 

= pA  — qA"^  (20) 

where  p * (ax2)  ^ (ax  - l)  and  q s b2/4a^x2 

In  a similar  manner,  4 can  be  determined,  and  its  average  found 
over  a cycle  of  The  result  is 

(«)av  - 0.  (21) 
If  the  average  value  of  A,  given  in  Eq.  (20),  is  used,  the 


variation  of  amplitude  with  time  can  be  found.  This  equation,  another 

example  of  Bernoulli* s equation,  can  be  made  linear  by  the  substitution 

—2 

w *=  A • The  resulting  linear  equation  is 


dw/dt  + 2pw  = -2q. 

(22) 

Its  solution 

is 

w ~ C exp(-2pt)  + q/p 

(23) 

and  thus, 

A « fc  exp(-2pt)  + q/p]“^2^ 

(24) 

where  C is 

an  arbitrary  constant.  If  A = Aq  at  t = 0, 

amplitude 

A is  given  as  a function  of  time  by 

A = A(t)  = jq/p  + (Aq  2 -*  q/p)  exp(~2ptfj 

(25) 

The  solution  just  obtained  is  only  approximate  because  just  the 
first  three  terms  of  the  infinite  series  of  Eq,  (19)  were  used.  This 


approximation  is  good  if  A/p  «1,  which  is  true  only  near  the  beginning 
of  the  growth  of  A and  provided  AQ<Cp,  As  the  instantaneous  amplitude 
increases,  the  solution  is  increasingly  in  error. 

Quantity  p in  Eq,  (25)  may  be  either  positive  or  negative,  as 
product  ax  is  greater  or  less  than  unity,  respectively.  If  p is 
positive,  amplitude  A approaches  a steady^state  value, 

= (f/q)1’'2  » 2p(ax  - i)^2 


(26) 


% 

16,  I 
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as  t becomes  infinite.  If  initially  Aq  is  small,  the  amplitude 
grows,  with  its  maximum  rate  of  change  occurring  when  X *=  0,  at  which 
instant  A «=  A Jj  , If  p is  negative,  amplitude  A vanishes  as 
t becomes  infinite. 

Plots  of  the  variation  of  A with  A,  and  of  A with  t,  as  1 

given  by  Eqs,  (20)  and  (25),  are  shown  in  Fig,  6,  Since  Eq,  (21) 

shows  that  the  average  frequency  is  constant,  the  approximate  solution 

for  Eq,  (9),  just  found,  is 

x = p + A(t)  sin(2^2t/T  + 9q)  (27) 

where  A(t)  is  given  by  Eq,  (25)  and  plotted  in  Fig,  6,  and  is  a 

phase  angle  determined  by  initial  conditions. 

Limitations  can  be  set  on  the  value  of  product  ax  allowed  fcr 

this  solution.  The  assumption  is  made  initially  that  the  change  in 

amplitude  per  cycle  of  the  oscillation  must  be  small.  If  p is 

positive,  the  maximum  value  A for  the  average  rate  of  change  of 

l/2 

amplitude  occurs  for  the  amplitude  Am  = (p/iq)^  , The  relative  rate 

of  growth  is  then,  from  Eq,  (20),  A ^ A^  = 2p/3»  The  period  of  the 

l/  2 

oscillation  is  T = 2rr/a  = 7T1  rrx.  Thus,  the  maximum  relative  change 

in  amplitude  per  cycle  is  nearly 

T Am/Am  - (21//2rrx)  (2p/3)  = (23//2tt/3)  (ax  - lKax)”1.  (2b) 

This  quantity  must  be  less  than  unity,  say,  tr  meet  the  assumption 

inherent  in  the  method  of  solution.  Thus,  the  requirement  is  that 

ax<[l  - (23^2n/3)"1]"1  - 3/2.  (29) 

If  p is  positive,  the  amplitude  approaches  the  steady  value, 

l/? 

Ag  * 23 (ax  - 1 )'v  This  amplitude  exceeds  the  value  3 for  ax>5/4, 
and  would  then  require  instantaneous  values  of  x to  change  algebraic 
sign.  Since  sign  changes  cannot  occur  because  of  infinite  damping  at 
x = 0,  a more  realistic  upper  limit  for  ax  in  the  solution  is  ax  <5/4. 


‘I 
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i_  p is  negative,  the  amplitude  decays  with  time.  According  to 
Eq.  (20),  the  maximum  rate  of  decay  occurs  for  the  largest,  or  initial, 
amplitude,  A rather  crude  estimate  of  the  lower  limit  for  ax  can  be 
obtained  by  disregarding  the  nonlinear  term  in  Eq,  (20),  in  which  case 
the  relative  race  of  decay  is  A/A  = p,  independent  of  the  amplitude. 
The  relative  change  in  amplitude  per  cycle  is,  then, 

T A/A  = (2^^rrx)(ax  - l)(ax2)  ^ . (30) 

This  quantity  must  be  greater  than  minus  one,  say,  so  that 

ax>  [l  + (21^2rr)'"1]'1  = 4/5.  (31) 

Thus,  in  order  that  the  approximate  solution  found  by  variation 
of  parameters  be  reasonably  accurate,  the  value  of  product  ax  must 
fall  within  the  limits 

4/5  < (ax)  <5/4. 

For  these  conditions,  the  approximate  solution  is  an  oscillation, 

l/2 

approximately  sinusoidal  with  angular  frequency  2 / /x,  having  its 
mean  value  at  x 53  |3.  Its  amplitude  decays  if  ax*'.  1.  The  amplitude 

1 /p 

approaches  a steady  value,  A = 2g(ax  - l)x/  , if  ax>l. 

s 

III.  4*  Approximate  solution  by  iteration,  about  (x  ,y  ) 

s s 

The  approximate  solution  obtained  by  the  method  of  variation  of 
parameters  consists  of  an  oscillation,  essentially  sinusoidal  in  wave- 
form, with  its  amplitude  varying  in  time.  An  oscillatory  steady  state 

is  achieved  if  ax  >1.  More  information  about  the  waveform  in  the 

13 

steady  state  can  be  found  by  a process  of  iteration. 

If  the  terms  in  x are  omitted  from  Eq.  (9),  it  becomes 

x + a2x  = a2[3  (32) 


13.  J.  J.  Stoker,  Nonlinear  Vibrations,  (Interscience  Publishers, 
New  York,  1950),  p,  83 
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which  has  the  solution 


x * (3  + A sin  at  (33) 

subject  to  the  initial  conditions  of  x * 3>  x « aA  at  t » 0,  This 
generating  solution  can  now  be  put  into  those  terms  of  Eq,  (9)  in 

which  x appears  so  as  to  give 

••2  2 3 

x + ax  = a p + a^xA  cos  at 

- (a3A/b3)  cos  at  (1  - a/3  sin  at 

+ A2/32  sin2  at  + • • - . ) (34) 

where  the  series  is  accurate  only  if  A/3^Cl«  Use  can  be  made  of 

the  identities 

sin  at  cos  at  « (-)  sin  2at 


2 1 

sin  at  cos  at  =■  (^)(cos  at  - cos  3at) 


to  give 

x + a2x  **  a2p  + a3A(x  - l/b3  - A2/4b33)  cos  at 

+ (a3A2/2b  32)  sin  2at  + • • • • (35) 

In  order  to  avoid  a secular  term  in  the  solution,  something  which  cannot 
occur  since  it  grows  indefinitely  with  time,  the  coefficient  of  cos  at 
must  be  zero,  or 

A - 23(ax  - l)^2.  (36) 

This  is  the  steady-state  amplitude  of  oscillation,  and  is  the  same 
result  found  in  Eq,  (26), 

Solution  for  Eq,  (35)  is,  then 

x e 3 + A sin  at  -»  aA2/6b32  sin  2at  + » * • * (37) 

where  the  coefficient  of  a possible  term  in  cos  at  is  made  zero  to 
keep  x = 3 at  t = 0, 

The  waveform  of  the  steady-state  oscillation  is  approximated  by 
Eq*  (37) j and  considerable  second-harmonic  distox’tion  is  seen  to  be 
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present.  The  relative  amplitude  of  the  second  harmonic  compared  with 
the  fundamental  component  is 

2nd/lst  =■  aA/6bp2  » (2^2/3)(eLX  - l)1/2^)"1.  (38) 

The  largest  value  of  ax  which  will  keep  instantaneous  x positive 
(as  is  necessary)  is  ax  =*  5/4*  For  this  value,  the  relative  second 
harmonic  is  2nd/lst  =■  1/5. 

A sketch  of  the  resulting  waveform  is  shown  in  Fig,  7.  The  presence 
of  the  harmonic  actually  makes  x go  negative  momentarily  in  this 
example.  The  rise  in  x from  values  just  greater  than  zero  to  large 
values  occurs  more  slowly  than  does  the  drop  from  large  values  to  small 
ones. 

It  should  be  recognized,  of  course,  that  Eq.  (37)  is  only  approxi- 
mate because  of  the  gross  assumptions  made  in  its  derivation.  It 
applies  reasonably  well  only  between  the  limits  l<ax<5/4, 

III,  5,  Solution  by  isocline  construction 

\ It  is  profitable  to  study  in  still  more  detail  the  solutions  of 

Eq,  (9)  as  represented  graphically  on  the  phase  plane.  In  carrying  the 

analysis  further  it  is  convenient  to  normalize  the  quantity  y “ x by 

defining  a new  variable  z h xy  *=  xx.  The  dimensions  of  z are  then 

the  same  as  those  of  x.  In  terms  of  z,  Eq,  (11)  becomes 

dz  = 2[(1  - l/bxx) z - x + g] 
dx  z 

This  equation  gives  the  slope  of  a solution  curve  at  any  point  \j  the 
z-x  plane.  If  this  slope  is  assigned  a constant  value,  say  m,  the 
locus  of  those  combinations  of  z and  x giving  this  assigned  slope 
can  be  found.  This  locus  is  the  isocline^  curve  connecting  points 

-n  ...  -h  ,,  „ ___  — — . - , _ 
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14.  N,  Minorsky,  ref.  7,  p.  20 


of  constant  slope,  m.  Its  equation  is 


- - 1 rstei  x)  . <«> 

Isocline  curves  for  any  assigned  value  of  m pass  through  the 

singular  point,  x = 8,  z = 0,  found  earlier.  All  curves  also 

s s 

pass  through  the  origin,  x 3 0,  z = 0,  which  is  another  singular 
point  of  more  complicated  nature  than  that  found  first.  The  slope  of 
all  isocline  curves  as  they  pass  through  the  origin  is  the  same, 
dz/dx  = ax, 

A set  of  isoclines  is  plotted  in  Fig,  8 upon  axes  of  z and  x» 

For  this  figure,  positive  numerical  values  of  the  parameters  are 

assumed  as  p=l,  bx  = 1,  so  that  also  ax  = 1,  According  to  the 

analysis  based  on  variation  of  parameters,  this  value  ax  = 1 is  just 

sufficient  to  lead  to  an  oscillating  solution  about  the  singular  point 

(x  , z ) with  the  amplitude  of  the  oscillation  remaining  constant, 
s s 

Isoclines,  calculated  from  Eq.  (40),  are  plotted  in  Fig.  8 for  several 
values  of  slope  m,  and  line  segments  drawn  through  the  isoclines 
have  the  corresponding  slope, 

Several  solution  curves  are  sketched  in,  always  cutting  the 
isoclines  with  the  required  slope.  The  predicted  oscillation  about 
the  singular  point  is  seen  to  occur.  Since  the  solution  curve  is  not 
a circle,  but  is  distorted,  the  solution  for  x vs.  t is  not  a 
simple  sinusoid,  although  it  is  periodic.  A solution  curve  coining 
from  a large  positive  value  of  x bends  so  as  to  approach  the  z-axis 
closely,  but  never  to  cross  it.  This  effect  comes  about  from  the  very 
large  positive  damping  occurring  when  x approaches  zero  from  the 
positive  side,  as  has  been  discussed  previously.  Because  of  this 
action,  a limit  cycle  appears.  For  the  assumed  values  of  bx  and  8, 


a positive  initial  value  of  x leads  to  an  oscillation  about 
the  singular  point. 

If  x is  initially  negative,  the  solution  curves  sketched  in 

Fig.  8 indicate  that  x always  returns  to  zero  with  a value  of  z 

that  goes  to  infinity.  Presumably  x will  then  pass  through  zero, 

with  x infinite  at  that  point,  and  become  positive.  As  soon  as  x 

is  positive,  the  solution  curve  spirals  inward  about  the  point  x = p, 

s 

z "0,  as  before.  It  appears,  therefore,  that  regardless  of  the 
s 

initial  algebraic  sign  of  x,  ultimately  it  becomes  positive. 

While  the  nature  of  the  singularity  at  the  origin  cannot  easily 

be  studied  so  completely  as  was  the  singularity  at  x =»  [3,  z =*  0, 

s s 

it  is  evident  in  this  example  that  the  origin  is  an  unstable  point. 

For  positive  values  of  x,  the  origin  is  similar  to  a saddle  point; 
for  negative  values  of  x it  is  similar  to  an  unstable  nodal  point. 

III.  6.  Approximate  solution  for  negative  x 

Some  of  the  qualitative  aspects  of  the  solution  obtained  for 
negative  values  of  x can  be  found  through  the  following  approximate 
analysis.  The  differential  equation  is 

(t^/2)  x - tx  + x/bx  + x = a/b  (41) 

which  is  merely  Eq.  (9)  with  the  coefficients  rearranged. 

If  t = 0,  this  equation  reduces  to  Eq.  (4),  having  the  exact 
solution  given  in  Eq.  (7),  This  latter  equation,  rewritten  here  for 
convenience  is 

X1  = |s  + h exp(-at'7j~1  (42) 

where  gab/a  and  h = (xq  - g),  with  x = xo  at  t = 0. 
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If  t >0,  but  is  small,  an  approximate  solution  for-  x can  be 
found  by  an  iteration  process.  The  solution  for  x1 , given  by  Eq.  (42) 
for  t = 0,  can  be  used  with  Eq,  (41)  for  x>0,  to  give  a new  solu- 
tion for  x as  follows.  From  Eq,  (41)  the  new  solution  is 

X£  = x^  + xx^  - (x2/2)  x^  (43) 

where  is  an  approximate  solution,  valid  for  x small.  In  writing 

this  equation,  is  put  in  place  of  the  third  and  fifth  terms  of 

Eq,  (41) , to  which  it  is  equivalent, 

A tentative  conclusion  from  the  isocline  plot  of  Fig,  8 was  that 
if  x initially  is  negative,  it  always  returns  to  aero,  even  if 
initially  x is  negative.  In  order  for  x to  return  to  zero,  x 
must  become  positive.  If  x is  initially  negative,  there  must  be 
some  instant  at  which  it  passes  through  zero  as  it  changes  sign. 

Thus,  the  time  at  which  x = 0 is  of  interest.  This  time  can  be 
determined  approximately,  making  use  of  Eq.  (43).  The  time  derivative 
of  this  equation,  with  the  assumption  that  x is  small  enough  that 
the  last  term  can  be  neglected,  is 

Xg  * + xx^.  (44) 

Derivatives  of  x., , found  from  Eq.  (42)  and  used  here,  give 

JIs-QL  - ax)  + h exp(-at)(l  + ax 
jg  + h exp(-at)j3 

If  X£  = 0,  the  requirement  is  that,  approximately, 

exp(-at)  = -(g/h)(l  - 2ax) 

= (l  - 2ax)/(l  - a/bxo).  (46) 

If  x = 0,  the  value  of  t from  this  equation  is  the  same  as  that 
found  earlier  as  the  time  for  which  x,  starting  with  a negative  value, 
goes  to  minus  infinity.  If  x is  made  slightly  greater  than  aero,  the 
time  for  0,  given  by  F,q,  (46),  is  also  made  larger. 


’ B ah  exp  (-at) 
*2 
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The  rate  at  which  x changes,  for  x ju3t  greater  than  zero  also 
can  be  estimated.  If  x = 0 and  if  the  magnitude  of  x is  small 
enough,  it  is  given  approximately  by 

x1  = xo  exp (at) . (47) 

This  equation  results  frcm  Eq.  (42)  if  jxoj<<,a/b.  This  value  of  x^ 
used  in  Eq,  (43)  leads  to 

(x^  - x1)/x1  = ax(l  - ax/2).  (48) 

This  equation  gives  the  difference  between  solution  x^,  applying 
for  x = 0,  and  solution  X£  applying  for  x >0,  but  small.  Since 
x^  and  x^  must  have  the  same  algebraic  sign,  the  left  side  of 
Eq.  (48)  is  positive  if  [x^l^Jx^j,  and  negative  if  JxgjcijjCjJ.  It 
is  evident  from  the  right  side  of  Eq.  (48)  that  the  fraction  is  positive, 


ard  thus  x, 
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These 


for  0<ax<rv2j  otherwise  jx,_ 
conclusions  apply  only  so  long  as  the  magnitude  of  x remains  small 
enough. 

The  result  of  thi3  discussion  can  be  summarized  in  Fig.  9,  in 

which  x is  plotted  against  t,  with  parameters  a and  b held 

constant.  The  initial  values  for  both  x = x and  x « x are  nega- 

o o 

tive,  and  both  x and  x are  assumed  the  same  for  all  curves.  The 
indication  is  that  unless  x = 0,  the  curves  for  x always  return  to 
zero,  and  that  the  point  where  x = 0 occurs  at  later  times  as  t is 
increased.  If  0<ax<2,  the  magnitude  jxj  initially  increases 
faster  than  for  t = 0.  If  ax >2,  the  magnitude  j x j initially 
increases  more  slowly  than  for  x **  0.  The  critical  value,  ax  = 2. 
separating  these  kinds  of  solutions  is  only  approximate  because  of  the 
crude  method  by  which  it  was  obtained. 


<xX  >2 


Fi^.S.  Qualitative  nature  of  solution  tor  Ec^.tS) 
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III.  7.  Solution  by  analog  computer 

The  differential  equation,  Eq.  (9),  was  set  up  on  a Reeves  Analog 
Computer,  using  the  system  diagrammed  in  Fig.  10.  Control  of  the 
parameters  a,  8 (*  a/b),  and  x is  possible,  as  shown.  The 
additional  inverting  amplifier  must  be  used  when  negative  values  of  x 
are  being  studied,  so  as  to  give  the  correct  algebraic  sign  to  the 
quotient  term.  The  computer  gives  accurate  results  so  long  as  x and 
its  derivatives  are  not  too  large.,  and  so  long  as  x is  not  too  small. 

If  x approaches  zero,  the  division  in  the  nonlinear  term  is  inaccurate. 

Some  typical  curves  representing  solutions  for  the  equation  are 

shown  in  the  accompanying  figures.  In  Figs.  11  and  12,  a plot  of  x 

against  t is  shown.  In  both  these  figures,  (3  = 20}  in  Fig,  11,  t 

is  varied  keeping  a constant,  while  in  Fig.  12,  a is  varied  keeping  t 

constant.  Initial  conditions  in  both  cases  are  xq  = +p/4  “ 5, 

xq  = yQ  “ °«  When  the  product  ax  is  unity,  an  oscillation  about 

x = p occurs,  with  its  amplitude  decaying  slowly.  The  approximate 

analysis  of  Sec.  III.  3.  predicts  a limit  cycle  for  this  value  of  ax, 

but  with  very  small  amplitude,  A limit  cycle  of  small  amplitude  may 

actually  occur}  it  is  difficult  to  decide  from  the  computer  solution. 

The  period  for  the  solution  with  ax  *=  1 is  quite  close  to  that 

1/2 

predicted,  T = 7T  ttx. 

If  product  ax  exceeds  unity,  the  limit  cycle  with  nonsinusoidal 
waveform  becomes  apparent.  The  waveform  is  similar  to  that  predicted 
in  Sec,  III.  4,  If  product  ax  is  small,  the  final  value  of  x is 
approached  monotonically-  Curves  of  x = x^  exp(at)  and  of  x,  as 
given  by  Eq,  (7)  for  x = 0,  are  plotted  in  Fig.  II, 
curves  as  x varies  between  zero  and  infinity. 


These  are  limiting 
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Captions  for  Figs,  11-15 

These  figures  are  all  solutions  for  Eq.  (9)  obtained  with  the 
analog  computer  set  up  of  Fig,  10,  The  conditions  for  the  figures 
are  as  follows, 

11,  a = 1/2,  3 = 20,  t as  indicated,  x^  = .5,  xo  = 0 

12,  x = 2,  3 e 20,  a adjusted  to  give  value  for  ax  as  indicated. 


13,  a = 1/2,  p = 20,  t adjusted  to  give  value  for  ax  as  indicated, 

x -set  to  value  indicated  by  circle,  x = 0 
o 3 o 

14*  a = l/lO,  3 = 10,  x as  indicated,  x = -5,  x = -0,75 

o o 

. a = 1/10,  3 - 10,  x as  indicated,  initial  conditions  as  indicated 

by  circle 


15 
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In  Fig.  13  are  shown  several  phase-plane  curves  for  y plotted 
against  x.  Cn  each  diagram  are  shown  solution  curves  for  several 
initial  values  of  xq,  always  with  yQ  = 0.  The  initial  points  are 
indicated  with  circles.  The  limit  cycle  for  at  = 1.4  is  evident, 
and  a limit  cycle  of  small  amplitude  may  occur  for  ax  * 1.  The  solu- 
tion curve  for  x = 0,  given  by  Eq.  (4),  is  plotted  for  comparison 
on  the  diagram  for  ax  = 0.35.  Here,  y jumps  quickly  from  its 
initial  zero  value  to  a value  near  that  which  would  have  to  exist 
with  the  case  of  x **  0. 

The  phase-plane  curve  obtained  here  for  ax  = 1 has  a different 
vertical  scale  from  that  of  Fig.  8,  since  the  vertical  coordinates 
are  z and  y,  respectively.  Otherwise  the  curves  are  very  similar. 

In  Fig.  14,  a plot  of  x against  t is  shown,  with  the  initial 
value  xQ  being  negative.  Parameters  were  chosen  as  a = l/l°, 

P “ 10,  so  that  b * 1/100*  The  initial  value  of  x is  xq  = -5. 

A curve,  calculated  from  Eq.  (?)  for  x = C is  shown  for  comparison. 
The  initial  slope  of  this  curve,  y = -0.75,  was  used  as  an  initial 
condition  in  obtaining  the  other  curves  also.  The  curve  for  x 51  5, 
ax  - 0,5,  is  seen  to  depart  from  zero  more  rapidly  than  the  curve  for 
x = 0}  the  curve  for  x = 20,  ax  *>  2,  leaves  zero  more  slowly. 
Phase- plane  curves  of  y plotted  against  x are  shown  in  Fig.  15, 
again  for  negative  values  of  x.  Once  more,  the  curve  for  x = 0, 
calculated  from  Eq.  (4),  is  shown  for  comparison.  If  x is  small, 
the  solution  curve  either  rises  or  falls  abruptly,  dependent  upon 
whether  the  initial  point  is  just  above  or  below  the  curve  for  x = 0, 
If  x is  larger,  a slower  change  occurs,  but  the  return  to  x = 0 is 
evident.  The  curves  of  Figs.  14  and  15  agree  with  the  qualitative 
conclusions  of  Sec.  111,6. 


' 4 
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IV.  Differential-Difference  Equation 

IV.  1,  Comparison  with  differential  equation 

The  original  nonlinear  differential-difference  equation,  3q,  (3)> 
was  replaced  by  a nonlinear  differential  equation,  Eq.  (9),  for  the 
analysis  of  Sec.  III.  The  differential  equation  wa3  obtained  by 
replacing  the  difference  term  in  the  original  equation  by  the  first 
three  terms  of  a Taylor's  series,  Eq.  (8).  The  two  equations  should 
yield  similar  solutions  so  long  as  the  Taylor's  series  approximates 
the  difference  term  well*  The  approximation  is  good  provided  the  time 
delay  t is  small  enough  that  the  terms  emitted  from  the  series  are 
small  compared  with  those  retained.  In  particular,  it  is  necessary 
that  (x^/6)  d^x/dt^  be  much  smaller  than  (t2/2)  d2x/dt^  , and 
so  on.  This  inequality  holds  if  x is  sufficiently  small.  If  t is 
larger,  it  is  unlikely  that  the  approximation  is  valid.  Solutions 
for  the  two  equations  then  would  be  expected  to  differ  by  significant 
amounts. 

An  estimate  of  the  magnitude  of  the  terms  can  be  found  from  the 
approximate  solutions  of  Eqs,  (36)  and  (37).  If  product  ax  is 
slightly  larger  than  unity,  approximately  x * p(l  + sin  at).  Then 
x = -pa2  sin  at  and  S?  = ~Pa^  cos  at  so  that  (r^/6)T  = 21//2p/3  and 

O 

(t*"/2)  x = P,  where  only  the  amplitudes  appear  in  the  last  two  rela- 
tions. Thus,  the  first  terra  omitted  in  the  series  has  an  amplitude 
about  one  half  that  of  the  last  term  retained,  and  the  approximation 
is  relatively  poor.  If  product  ax  is  larger,  the  solution  changes 
more  abruptly  with  time,  and  higher-order  derivatives  are  larger. 

The  approximation  could  be  expected  to  be  even  less  accurate. 
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Furthermore,  the  differential-difference  equation  is  equivalent 
to  a differential  equation  of  infinite  order,  which  may  have  an  infinity 
of  modes  of  oscillation.  The  approximate  differential  equation,  of 
second  order,  can  have  only  a single  mode  of  oscillation.  It  is  not 
difficult  to  see  how  the  various  modes  for  the  difference  equation 
can  arise. 

The  analysis  of  the  differential  equation  has  predicted  that  if 

product  ax  exceeds  unity,  a limit  cycle  appears,  leading  to  an 

oscillation  of  period  T as  shown  in  Fig,  16.  This  solution  can  be 

assumed  to  apply  to  the  differential-difference  equation,  also. 

According  to  this  latter  equation,  the  value  of  x at  some  tine  t 

depends,  in  part,  upon  its  value  at  the  earlier  time  (t  - x) . If 

for  example,  the  time  (t  - x)  is  zero,  time  t at  which  x is 

determined,  can  be  t = x,  (x  + T),  (x  + 2T),  * • • • (x  + nT) , 

For  any  of  these  times,  the  same  result  must  be  obtained  for  x, 

because  of  the  periodicity  of  the  solution.  The  same  steady-state 

solution  results  if  the  original  delay  time  is  increased  by  any 

integral  multiple  of  the  period.  The  solution  of  Fig,  16,  then,  can 

be  obtained  for  fixed  values  of  parameters  a and  b,  and  an  infinity 

of  values  for  x.  The  period  and  waveform  remained  unchanged. 

If  all  three  parameters,  a,  b,  and  x,  are  fixed,  a variety  of 

steady-state  solutions  also  may  occur.  An  example  is  shown  in  Fig.  17, 

where  the  waveforms  are  merely  sketched  roughly.  If  the  product  ax 

is  large  enough,  a violent  oscillation  of  long  period,  T^  = x/k,  may 

occur,  where  k is  a constant.  If  ax  is  just  larger  than  unity, 

1/2 

the  period  of  the  oscillation  has  been  shown  to  be  T - 2 ttt,  so 
that  k * (2^^rr)  1 = (4,44)  \ Fcr  larger  values  of  ax,  k vail 
be  smaller,  but  the  effective  time  delay  remains  in  the  order  of  l/4 
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of  a period.  For  an  initial  value  of  ax  sufficiently  large,  another 
steady-state  oscillation  may  occur,  with  period  =*  x/(l  + k).  The 
effective  time  delay  is  then  in  the  order  of  5/4  a period,  and  the 
period  of  the  second  oscillation  is  about  l/5  that  of  the  first.  The 
waveform  of  the  second  oscillation  is  determined  by  the  product  ax* , 
where  the  value  of  x<  is  the  least  delay  that  could  give  the 
observed  oscillation,  or  about  l/5  the  actual  value  of  x.  Thus,  the 
waveform  of  the  oscillation  of  shorter  period  depends  upon  the  quantity 
ax/5,  approximately,  and  if  there  is  to  be  a steady-state  oscillation 
this  quantity  must  be  large  enough  to  lead  to  a limit  cycle,  that  is, 
somewhat  larger  than  unity.  The  oscillation  of  longer  period  vdll  have 
its  waveform  dependent  upon  ax,  which  will  be  greater  than  five,  and 
the  waveform  will  involve  large  peaks  and  violent  changes. 

If  product  ax  were  sufficiently  large,  a whole  sequence  of 
steady-state  oscillations  of  this  sort  might  be  obtained,  each  having 
a different  period  and  a different  waveform.  Suitable  and  rather 
special  initial  conditions  would  be  required  to  start  3uch  oscillations, 
however.  It  is  likely  that  only  the  oscillation  of  longest  period 
will  actually  occur  in  any  physical  system  to  which  the  equation 
applies, 

IV.  2,  Solution  by  analog  computer 

The  differential-difference  equation  is  difficult  to  set  up  on  a 
simple  analog  computer  because  of  the  necessity  for  introdixeing  the 
time  delay  x.  This  delay  time  would  have  to  be  reasonably  long  in 
comparison  with  the  period  of  oscillation  of  the  highest  frequency 
that  can  be  handled  by  the  computer.  It  is  not  easy  to  set  up  a system 
with  the  available  components  so  as  to  obtain  the  necessary  delay. 
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In  order  to  avoid  the  need  of  introducing  a time  delay  in  the 
computer  itself,  a step-by*-step  computation  was  used.  The  computei 
was  used  to  solve  the  equation 

dx(t)/dt  = jji  - b x(tQ  + At/2  - v)]  x(t)  (49) 

in  place  of  Eq.  (3).  In  Eq.  (49),  At  is  the  time  interval  during 
which  computation  occurs  and  t is  the  tine  at  the  beginning  of  the 
interval.  In  the  bracket,  a constant  value  for  x is  used,  evaluated 
at  the  delay  time  r earlier  than  the  average  of  the  time  during  the 
interval.  If  the  time  interval  is  chosen  small  enough,  this  procedure 
should  lead  to  accurate  results.  It  suffers,  of  course,  from  the 
usual  errors  of  step~by*-step  computations. 

The  computer  was  set  up  as  shown  in  Fig.  18,  The  procedure  is  to 
put  into  the  multiplying  circuit  of  the  computer  the  value  of 
x(tQ  + At/2  "•  t).  The  value  of  x(tQ)  is  put  into  the  integrator  as 
an  initial  condition.  The  computer  is  allowed  to  run  for  the  time 
interval  At,  after  which  the  value  of  x(tn  + At)  is  read.  This 
new  value  of  x is  then  used  for  new  initial  conditions  with  the 
next  time  interval,  and  the  procedure  is  repeated.  In  this  way 
successive  points  for  a curve  of  x versus  t can  be  obtained. 

A family  of  curves  found  in  this  way  is  shown  in  Fig,  19,  where 
x is  always  positive.  Numerical  values  for  the  parameters  were  t = 18, 
8 = 20,  and  At  = 4.  Parameter  a was  adjusted  to  give  product  av 
the  indicated  value.  The  initial  value  of  x at  t = 0 was  taken  as 
xq  * 2.  It  was  assumed  further  than  x had  this  same  value  during 
the  four  intervals  At  preceding  t = 0.  This  assumption  gives  x 
the  constant  value,  x = 2,  during  the  time  -16  < t < 0,  after 
which  x begins  to  change,  With  a value  of  x as  small  as  this. 
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Fi^.  18.  Analog  computer  setup  for  solving  Etj.  14S) 
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Captions  for  Figs.  19-20 

These  figures  are  solutions  for  Eq.  (3)  obtained  vd.th  the  analog 
conputer  set  up  of  Fig,  18.  The  conditions  for  the  figures  are  as 
fallows. 

19.  t = 18,  0 * 20,  At  3;  4,  a adjusted  to  give  value  for  at  as 

indicated,  x = 2 for  -16  t < 0 

20,  t ■ 18,  0 = 20,  At  = 4,  a adjusted  to  give  value  for  ax  as 

indicated,  x * -2  for  -16  t < 0 


aT*2.0 


in  comparison  with  3 = 20,  the  exact  value  used  in  starting  the 
computation  does  not  influence  the  resulting  solution  curve  very  much. 

The  curves  of  Fig,  19  for  the  differential-difference  equation 
are  similar  qualitatively  to  those  of  Fig.  12  for  the  differential 
equation.  The  period  of  the  oscillation  is  similar  for  a given  value 
of  product  ax.  However,  the  minimum  value  of  ax  leading  to  a 
limit  cycle  for  the  differential  equation,  ax  » 1,  is  less  than  that, 
ax  * 1.6,  needed  for  the  differential-difference  equation,  A differ- 
ence of  this  sort  is  not  unexpected,  however,  because  of  the  approxi- 
mation that  was  used  in  getting  the  differential  equation.  In  the 
appendix  i3  given  a discus  S 1 On  VrlijL  ch  might  indicate  the  limit  cycle 
would  first  appear  near  ar  = tt/2. 

Another  family  of  curves,  with  x negative,  is  shown  in  Fig,  20. 
The  same  numerical  values  of  the  parameters  were  used  in  obtaining 
Fig.  20  as  in  Fig.  19,  except  for  the  initial  value  of  x.  For  Fig,  20 
it  was  assumed  that  x = -2  for  the  time  -16  < t < 0.  The  curves  of 
Fig.  20  show  the  expected  shape,  with  x becoming  an  increasingly 
large  negative  quantity  as  t increases.  The  rate  of  increase 
depends  upon  the  value  of  parameter  a.  The  monotonic  increase  of 
the  curves  of  Fig.  20  is  again  different  from  the  solutions  for  the 
differential  equation,  shown  in  Fig,  14,  where  x always  returns  to 
zero.  Again,  this  results  because  of  a basic  difference  in  the  two 
equations. 

V.  Conclusion 

The  differential-difference  equation  has  been  approached  in  a 
variety  of  ways,  and  a number  of  approximate  solutions  for  it  have  been 
obtained.  Most  complete  information  comes  frcm  the  differential  equa- 
tion that  is  approximately  equivalent.  Information  about  the  solution 
is  summarized  in  Tables  1 and  2. 
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Table  1 

Differential-Difference  Equation 

1.  Form:  dx(t)/dt  = ja  - b x(t  - t)"|  x(t)  (3) 

a,  b,  x are  positive  real  constants 

2.  Equilibrium  conditions: 

x(t)  **  x(t  — t)  =0,  unstable 

x(t)  = x(t  — t)  *»  a/b,  stable  with,  perhaps,  a limit  cycle 

3.  Algebraic  sign  of  solution:  solution  cannot  change  sign 

4..  Exact  solution  for  x = 0: 

x = [b/a  + (x^1  - b/a)  exp(-a  (7) 

5.  Approximate  solution  for  x>0:  qualitatively  similar  to  solutions 

for  the  differential  equation,  Eq.  (9).  See  Fig,  19. 

6.  Probable  value  of  ax  for  oscillation,  x>0:  ax>  l/2, 
this  value  is  estimated 

7.  Probable  value  of  ax  for  a limit  cycle,  x>  0:  ax>n/2 

a.  Solution  for  x<0;  solution  always  gees  to  x = ~co  , See  Fig,  20. 
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Differential  Equation 

1.  Form:  5c  - a,2 x*  + a2  %/hx.  + a2  x = a2p  (9) 

2 2 

a = 2 :/x  , 3 = a/b,  where  a,  b,  t are  positive  real  constants 

2.  Equilibrium  conditions: 

x a 0,  x = 0,  unstable 

x = a/b,  x = 0,  stable  with,  perhaps,  a limit  cycle 

3.  Algebraic  sign  of  solution:  Starting  with  either  sign,  solution 

always  becomes  positive. 

4*  Exact  solution  for  x = 0: 

x = fb/a  + (xq  1 - b/a)  expC-atyi"1  (?) 

5.  Approximate  solution  for  x>0: 

For  4/5<ax<5/4: 

x - p + A(t)  sinte^t/x  + eQ)  (2:') 

A(t)  = [q/p  + (Aq~2  - q/p)  axp(-2pt)j  1^2  (25) 

p = (ax2)  (ax  - 1),  q = b2/4a^x2 

at  t = 0,  A * A , $ = ■©  . 

’ o’  o 

In  steady  state,  l<ax<C5/4: 

x = p + A sin  at  - aA  /6bp  sin  2at  + • • • • (37) 

See  Fig,  12 

6.  Probable  value  of  ax  for  oscillation,  x>0: 

ax>(2^2  - 1)  = 0,414 

7.  Probable  value  of  ax  for  a limit  cycle,  x>0: 

ax  > 1 

8.  Solution  for  x<0:  solution  always  ultimately  returns  to 

x = 0 and  becomes  positive,  with  x = uo  at  x = 0,  Nature  of 
solution  depends  upon  parameters  of  equation  and  initial  conditions. 
See  Fig.  14, 
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Appendix  A.  Method  of  Equivalent  Linearization 

There  is  yet  another  method  that  is  sometimes  useful  in  finding 

an  apprcocimate  solution  for  an  equation  that  cannot  be  solved  exactly. 

This  is  the  method  of  equivalent  linearization,  sometimes  referred  to 

lr  16 

under  the  names  of  Kryloff  and  Bogoliuboff.  It  can  be  applied 

to  the  equation  under  discussion  here,  wi  ch  is 

dx(t)/dt  s j^a  - b x(t  - t)  x(t).  (3) 

Under  some  conditions  this  equation  has  been  shown  to  have  steadjr- 
state  solutions,  oscillating  about  a mean  value. 

The  procedure  is  to  assume  a solution  of  the  necessary  form 
I x(t)  = P + Q sin  cot  (A.l) 

where  P,  Q,  and  u)  are  all  constants  to  be  determined.  The 
• assumed  solution  is  then  substituted  into  Eq.  (?),  giving 

t 

t 2 i o 

Qco  cos  cot  a (aF  - bP  - ~ bQ~  cos  rot) 

+ (aQ  - bFQ  - bFQ  cos  cot)  sin  cot 

+ (bPQ  sin  cox)  cos  cot 
1 2 

+ 2 hQ  (cos  WT  003  2cot  + sin  cot  sin  2cot)« 

It  is  argued  that  this  equation  must  be  valid  for  those  components  of 
zero  frequency  and  the  fundamental  frequency,  anu  terms  of  frequency 
higher  than  the  fundamental  are  ignored.  The  following  three  relations 
are  thus  obtained. 

2 12  ■\ 
constant:  aP  - bP  - — bQ  cos  cot  = 0 > 

i 

sin  cot:  aQ  - bPQ  - bPQ  cos  cot  - 0 j (A. 2) 

cos  cot:  Qjj  - bPQ  sin  cot  = 0 ) 

These  equations  are  not  correct  if  t = 0,  since  then  co  = 0,  and 
the  term  in  cos  2cot  becomes  a constant  and  must  be  considered. 

15 « N*  Mi'.icrsky,  ref,  7,  Ch,  XII 

<* . 

16.  F,,  E.  Bothwell,  Econometrica,  20,  269,  (1952) 
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The  relations  can  be  rearranged  to  give 

£ 2 

sin  (jjt  = 2ato/(a  + ur) 

cos  cut  * (a2  - u:2)/(a2  + <u2) 

P = (a2  + o)2)/2ab 

Q - 2^2  P.  (A, 5) 

These  equations  must  be  satisfied  simultaneously,  and  can  be  solved 

for  the  necessary  values  of  o,  P,  and  Q» 

It  is  evident  from  the  transcendental  form  of  Eq.  (A, 3)  that  an 

infinity  of  values  of  u>,  the  angular  frequency  of  oscillation,  are 

allowed.  This  is  the  result  expected,  since  the  original  equation  is 

equivalent  to  a differential  equation  of  infinite  order.  However,  the 

relation  between  P and  Q,  Eq.  (A. 5),  cannot  correspond  to  a correct 

solution  for  Eq.  (3).  It  requires  that  Q exceed  P,  so  that 

instantaneous  values  of  x always  change  sign  twice  during  a cycle  of 

the  oscillation,  something  which  cannot  occur.  Furthermore,  Eq.  (A. A) 

indicates  that  the  amplitude  of  oscillation  is  larger  for  a higher- 

frequency  mode  of  oscillation,  which  seems  to  be  unreasonable r 

If  product  ax  = tt/2,  Eq.  (A. 3)  predicts  that  = n/2  13  1.57, 

and  Eq.  (A. 4)  predicts  that  P-^  * a/b.  This  result  agrees  fairly 

well  with  the  discussion  of  Sec.  III.  3 where  for  ax  near  unity,  it 

1/2 

was  found  that  wx  = 2 K 1.41  and  the  mean  value  of  the  solution 
is  0 = a/b.  A second  value  of  wx,  allowed  by  Eq.  (A.3)  with 
ax  = n/2,  is  about  w2x  = 9*1  or  approximately,  w2  = 5.S  This 
result  Is  also  about  what  is  expected,  with  the  ratio  of  the  first  two 
frequencies  of  oscillation  being  in  the  order  of  five  to  one.  The 
value  of  P for  this  second  frequency,  from  Eq.  (A. 4),  is  predicted  as 
P2  = 17  a/b,  and  this  seems  far  too  large.  Furthermore,  numerical 


(A.3) 


ft.  1 \ 
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value s found  from  Eas.  (A -.3**5 ) do  not  appear  reasonable  for  values  of 
product  ax  differing  from  ax  = tt/2. 

The  conclusion  must  be  drawn  that  this  method  of  solution  yields 
results  of  seme  utility  only  for  a narrow  range  of  values  of  product 
ax  near  ar  « tt/2.  Most  likely  this  is  about  the  value  of  product  ax 
which  first  gives  a limit  cycle  in  the  solution  of  the  differential- 
difference  equation.  In  general,  however,  the  method  gives  results  that 
seem  quite  unreasonable.  Presumably  this  occurs  because  solutions 
for  Eq.  (3)  are  sufficiently  non-sinusoidal  that  the  assumed  form  of 
Eq.  (A.l)  is  far  from  correct. 


r 


Appendix  B.  Aj 


Solution  for  a~r  !»1 


An  exact  solution  for  Eq.  (3)  with  the  time  delay  equal  to  zero 
was  found  in  Sec.  II.  Solutions  with  values  of  product  ax  in  the 
order  of  unity  were  found  in  Sec.  IV.  It  is  of  interest  to  explore 
the  nature  of  the  solution  if  ax  becomes  very  large. 

It  was  shown  in  Secs.  Ill  and  IV  that  if  ar  is  sufficiently 
large,  the  solution  is  periodic  with  the  period  T.  The  oscillations 
are  nearly  sinusoidal  for  ax  just  large  enough  to  give  the  periodic 
solution,  and  become  violent  relaxation  oscillations  for  av  very 
large.  Always,  however,  the  mean  value,  X,  of  the  solution  is  the 
same,  X » a/b  * 0.  This  relation  can  be  shov/n  as  follows. 


Equation  (3)  can  be  vritten  as 


tato/dt 


a - b x(t  - t). 


(B.l) 


If  the  solution  is  periodic  with  period  T,  both  sides  of  this  equa- 


tion can  be  averaged  over  a period  to  give  the  relation 

T T 

i / dx(t)  i fr  . n „ 

T I T La  - b - T'Jut> 


(B.2) 


Because  of  the  periodicity,  the  integral  on  the  left  side  of  the 
equation  is  zero.  Also,  the  mean  value  of  the  oscillating  solution  is 


X * (1/T)  / x(t)  dt  « (l/l)  ( x(t  - t)  dt. 


(B.3) 


Thus,  Eq.  (B.2)  becomes 


0 « a - b x 


This  fact  was  observed  by  L.  Onsager  (ref.  3,  P.  237).  The  deriva- 
tion given  here  was  suggested  by  P,  M.  Schultheies, 
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X = a/b  - p.  (BJ:) 

This  equation  is  valid  for  any  periodic  solution. 

As  product  at  is  made  large,  a plot  of  the  solution  for  Eq.  (3) 

approximates  the  shape  sketched  in  Fig.  B.l,  v/here  the  time  origin  is 

choEien  at  the  point  at  which  x crosses  the  value  x = p as  it 

increases  in  magnitude.  The  rise  from  x = p to  the  maximum  value, 

x , takes  place  almost  linearly.  The  drop  from  x „ to  a very 
max  max 

small  value,  nearly  x = 0,  takes  place  very  abruptly.  The  solution 
stays  at  a small  value  most  of  the  remainder  of  the  period,  and  then 
begins  to  rise  again  rather  suddenly.  The  shape  of  Fig.  B.l  is  consid- 
erably simplified,  of  course*  the  actual  solution  curve  does  not  change 
value  dis continuously. 

At  the  time  t = 0 in  Fig.  B.l,  x = p,  and  the  value  of  x at 

the  earlier  time,  (t  - x)  =»  (-x),  is  essentially  zero.  Thus,  at 

t «■  0,  the  slope  of  the  solution  curve  is  almost  dx(0)/dt  = a|3,  from 

Eq.  (3)*  It  is  assumed  that  this  slope  is  maintained  over  the  inverval 

0 K t <1  x.  At  the  end  of  the  interval,  t = x,  and  from  Eq.  (3), 

the  slope  must  vanish,  dx(x)/dt  » 0,  since  x(t  - x)  =■  3,  Thus,  at 

t * x,  x has  its  maximum  value,  x « kp,  this  relation  defining 

' max 

the  constant  k.  The  combination  of  the  known  slope  for  the  solution 
curve,  maintained  over  the  known  time  interval,  allows  the  value  of  k 
to  be  determined  as 

dx/dt  = ap  = (k  - l)p/x 


or 

k = (1  + ax). 

Thu3,  the  maximum  value  of  x is 


x 

vr 

Jl 


’.ax 


(B.5) 


(E.6) 


(1  + ax)p 


Since  it  was  shown  in  Eq,  (B.4)  that  the  mean  value  of  the 
periodic  solution  is  constant,  the  positive  and  negative  areas  between 
the  solution  curve  and  x = x “ 3 must  be  the  same  in  magnitude* 

Thus,  the  relation  can  be  written 

x(k  - 1)3/2  = (T  - t)3 

or  the  period  is 

T = (k  + l)x/2  = (1  + aT/2)x.  (B.7) 

Equations  (B.6)  and  (B,7)  c an  be  used  for  sketching  the  approximate 
shape  of  the  solution  curve,  provided  product  ax  is  known  and  is 
3’officiently  large.  Probably  the  maximum  value  of  x,  given  by 
Eq,  (B.6),  is  reasonably  accurate  for  ax  > 2.  The  period,  given  by 
Eq,  (B.7),  probably  is  not  very  accurate  unless  ax  > 10,  however. 

The  poor  accuracy  in  this  latter  equation  results  from  the  assumption 
of  discontinuous  jump?  in  the  solution,  and  these  approximate  the 
actual  solution  well  only  if  ax  is  quite  large. 

An  estimate  of  che  ratios  (x max/3)  and  (T/x)  for  periodic 
solutions  of  Eq,  (3)  is  shewn  in  Fig,  B.2.  These  estimates  are  based 
on  the  solutions  found  in  Secs,  III  and  IV  for  ax  near  two,  and  upon 
the  equations  just  developed  for  ax  very  large.  Because  of  the 
wide  range  in  values  of  x which  occurs  when  ax  is  lai’ge,  it  has 
not  been  possible  to  check  these  estimates  quantitatively  upon  the 
analog  computer. 
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Fig.  B.l  Approximate  solution  for  Eq.  (3)  with  ax  very  large. 
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Fig.  B.2  Estimated  period  and  maximum  value  of  solution  for  Eq.  (3). 


